Maximizing as minimum
Adapted from: (Floudas et al., 1999; Section 4.10), (Laurent, 2008; Example 6.23) and (Lasserre, 2009; Table 5.1)
Introduction
Consider the polynomial optimization problem from (Floudas et al., 1999; Section 4.10) of minimizing the linear function $-x_1 - x_2$ over the basic semialgebraic set defined by the inequalities $x_2 \le 2x_1^4 - 8x_1^3 + 8x_1^2 + 2$, $x_2 \le 4x_1^4 - 32x_1^3 + 88x_1^2 - 96x_1 + 36$ and the box constraints $0 \le x_1 \le 3$ and $0 \le x_2 \le 4$,
using DynamicPolynomials
@polyvar x[1:2]
p = -sum(x)
using SumOfSquares
f1 = 2x[1]^4 - 8x[1]^3 + 8x[1]^2 + 2
f2 = 4x[1]^4 - 32x[1]^3 + 88x[1]^2 - 96x[1] + 36
K = @set x[1] >= 0 && x[1] <= 3 && x[2] >= 0 && x[2] <= 4 && x[2] <= f1 && x[2] <= f2Basic semialgebraic Set defined by no equality
6 inequalities
x[1] ≥ 0
3 - x[1] ≥ 0
x[2] ≥ 0
4 - x[2] ≥ 0
2 - x[2] + 8*x[1]^2 - 8*x[1]^3 + 2*x[1]^4 ≥ 0
36 - x[2] - 96*x[1] + 88*x[1]^2 - 32*x[1]^3 + 4*x[1]^4 ≥ 0
As we can observe below, the bounds on x[2] could be dropped and optimization problem is equivalent to the maximization of min(f1, f2) between 0 and 3.
xs = range(0, stop = 3, length = 100)
using Plots
plot(xs, f1.(xs), label = "f1")
plot!(xs, f2.(xs), label = "f2")
plot!(xs, 4 * ones(length(xs)), label = nothing)
We will now see how to find the optimal solution using Sum of Squares Programming. We first need to pick an SDP solver, see here for a list of the available choices.
import Clarabel
solver = Clarabel.OptimizerClarabel.MOIwrapper.OptimizerA Sum-of-Squares certificate that $p \ge \alpha$ over the domain S, ensures that $\alpha$ is a lower bound to the polynomial optimization problem. The following function searches for the largest lower bound and finds zero using the dth level of the hierarchy`.
function solve(d)
model = SOSModel(solver)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c, p >= α, domain = K, maxdegree = d)
optimize!(model)
println(solution_summary(model))
return model
endsolve (generic function with 1 method)The first level of the hierarchy gives a lower bound of -7`
model4 = solve(4)-------------------------------------------------------------
Clarabel.jl v0.11.1 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 31
constraints = 40
nnz(P) = 0
nnz(A) = 73
cones (total) = 6
: Zero = 1, numel = 10
: PSDTriangle = 5, numel = (6,6,6,6,6)
settings:
linear algebra: direct / qdldl, precision: 64 bit (1 thread)
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 1.2590e-01 1.2590e-01 8.33e-17 7.72e-01 5.10e-01 1.00e+00 1.51e+00 ------
1 4.8332e-01 5.0333e-01 2.00e-02 1.61e-01 9.82e-02 1.86e-01 3.32e-01 8.02e-01
2 1.3995e+00 1.5138e+00 8.17e-02 8.80e-02 3.22e-02 2.10e-01 1.08e-01 7.85e-01
3 2.7854e+00 2.9400e+00 5.55e-02 1.81e-02 4.48e-03 1.83e-01 1.58e-02 8.78e-01
4 4.0242e+00 4.2592e+00 5.84e-02 1.24e-02 1.85e-03 2.61e-01 7.20e-03 7.46e-01
5 5.4090e+00 5.5096e+00 1.86e-02 3.46e-03 4.31e-04 1.09e-01 1.81e-03 8.27e-01
6 6.2811e+00 6.3385e+00 9.13e-03 1.45e-03 1.50e-04 6.15e-02 6.79e-04 7.22e-01
7 6.9685e+00 6.9720e+00 5.12e-04 5.52e-05 5.40e-06 3.75e-03 2.52e-05 9.69e-01
8 6.9993e+00 6.9994e+00 1.29e-05 1.39e-06 1.35e-07 9.46e-05 6.33e-07 9.75e-01
9 7.0000e+00 7.0000e+00 3.48e-07 3.79e-08 3.68e-09 2.56e-06 1.72e-08 9.73e-01
10 7.0000e+00 7.0000e+00 1.18e-08 1.30e-09 1.25e-10 8.65e-08 5.85e-10 9.66e-01
11 7.0000e+00 7.0000e+00 2.34e-10 2.62e-11 2.59e-12 1.72e-09 1.17e-11 9.80e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 1.63ms
solution_summary(; result = 1, verbose = false)
├ solver_name : Clarabel
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 1
│ └ raw_status : SOLVED
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ ├ objective_value : -7.00000e+00
│ └ dual_objective_value : -7.00000e+00
└ Work counters
├ solve_time (sec) : 1.63231e-03
└ barrier_iterations : 11The moment matrix is not flat so we cannot extract a minimizer with atomic_measure. We can still look for a feasible solution close to the vector of first-order moments $(\mathbb{E}[x_1], \mathbb{E}[x_2])$, as suggested in (Laurent, 2008). This vector is not necessarily feasible so round_solution projects it onto K with heuristic_projection. The constraints 0 ≤ x[1] ≤ 3 only depend on x[1] so x[1] is first projected onto [0, 3]. Once x[1] is fixed, the remaining constraints only depend on x[2] so x[2] is then projected onto [0, min(4, f1(x[1]), f2(x[1]))].
ν4 = moment_matrix(model4[:c])
x4 = round_solution(ν4, K, p)2-element Vector{Float64}:
2.9999999972673947
0.0The objective value at this feasible point is an upper bound to the optimal objective value so we now know that it is between -7 and -3.
p(x4)-2.9999999972673947Instead of only considering the first-order moments, we can sample from the Gaussian distribution with the same moments of order up to 2, project the samples and keep the best feasible point found. This generalizes the random hyperplane rounding of (Goemans and Williamson, 1995), see also (Barak and Steurer, 2016). We fix the seed of the random number generator so that the results are reproducible.
import Random
gaussian = GaussianRounding(rng = Random.MersenneTwister(0))
p(round_solution(ν4, K, p, rounding = gaussian))-3.0The second level improves the lower bound
model5 = solve(5)-------------------------------------------------------------
Clarabel.jl v0.11.1 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 108
constraints = 128
nnz(P) = 0
nnz(A) = 266
cones (total) = 7
: Zero = 1, numel = 21
: Nonnegative = 1, numel = 2
: PSDTriangle = 5, numel = (21,21,21,21,21)
settings:
linear algebra: direct / qdldl, precision: 64 bit (1 thread)
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 2.8313e-01 2.8313e-01 0.00e+00 8.12e-01 7.55e-01 1.00e+00 1.82e+00 ------
1 5.4905e-01 5.6026e-01 1.12e-02 1.47e-01 1.42e-01 1.86e-01 4.16e-01 8.05e-01
2 1.0075e+00 1.0482e+00 4.04e-02 3.92e-02 4.05e-02 8.87e-02 1.08e-01 8.48e-01
3 1.7352e+00 1.7929e+00 3.33e-02 8.93e-03 6.29e-03 7.10e-02 1.55e-02 9.20e-01
4 2.2703e+00 2.3425e+00 3.18e-02 3.95e-03 1.67e-03 7.97e-02 4.20e-03 8.20e-01
5 2.5204e+00 2.5826e+00 2.47e-02 2.72e-03 9.42e-04 6.81e-02 2.42e-03 6.21e-01
6 3.3688e+00 3.4392e+00 2.09e-02 7.06e-04 1.40e-04 7.25e-02 3.73e-04 8.78e-01
7 4.1613e+00 4.2336e+00 1.74e-02 3.97e-04 3.95e-05 7.36e-02 1.14e-04 8.80e-01
8 4.9013e+00 4.9482e+00 9.56e-03 1.79e-04 1.33e-05 4.76e-02 4.21e-05 8.78e-01
9 5.7265e+00 5.7548e+00 4.95e-03 6.71e-05 3.22e-06 2.86e-02 1.12e-05 8.50e-01
10 5.7843e+00 5.8135e+00 5.04e-03 6.48e-05 2.90e-06 2.95e-02 1.03e-05 2.28e-01
11 6.4423e+00 6.4536e+00 1.75e-03 1.41e-05 5.51e-07 1.13e-02 2.07e-06 8.34e-01
12 6.6492e+00 6.6500e+00 1.19e-04 8.39e-07 3.36e-08 7.94e-04 1.26e-07 9.69e-01
13 6.6662e+00 6.6663e+00 3.10e-06 2.18e-08 8.74e-10 2.08e-05 3.28e-09 9.74e-01
14 6.6667e+00 6.6667e+00 8.84e-08 6.27e-10 2.49e-11 5.93e-07 9.36e-11 9.72e-01
15 6.6667e+00 6.6667e+00 2.63e-09 1.88e-11 7.39e-13 1.76e-08 2.79e-12 9.71e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 5.70ms
solution_summary(; result = 1, verbose = false)
├ solver_name : Clarabel
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 1
│ └ raw_status : SOLVED
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ ├ objective_value : -6.66667e+00
│ └ dual_objective_value : -6.66667e+00
└ Work counters
├ solve_time (sec) : 5.69689e-03
└ barrier_iterations : 15The upper bound obtained from the first-order moments is improved as well:
ν5 = moment_matrix(model5[:c])
x5 = round_solution(ν5, K, p)
p(x5)-3.9012348317040884With the Gaussian rounding, we obtain a better upper bound:
x5_gaussian = round_solution(ν5, K, p, rounding = gaussian)
p(x5_gaussian)-5.424978206364312The third level finds the optimal objective value as lower bound...
model7 = solve(7)-------------------------------------------------------------
Clarabel.jl v0.11.1 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 288
constraints = 323
nnz(P) = 0
nnz(A) = 739
cones (total) = 8
: Zero = 1, numel = 36
: PSDTriangle = 7, numel = (55,55,55,55,...,55)
settings:
linear algebra: direct / qdldl, precision: 64 bit (1 thread)
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 3.3234e-01 3.3234e-01 3.89e-16 8.45e-01 6.51e-01 1.00e+00 1.51e+00 ------
1 4.7370e-01 4.8104e-01 7.34e-03 1.45e-01 1.49e-01 1.77e-01 3.55e-01 8.03e-01
2 7.3313e-01 7.6805e-01 3.49e-02 4.83e-02 6.81e-02 1.12e-01 1.54e-01 7.48e-01
3 1.1439e+00 1.1684e+00 2.14e-02 8.88e-03 1.34e-02 4.09e-02 2.71e-02 8.65e-01
4 1.5684e+00 1.6151e+00 2.98e-02 3.54e-03 4.11e-03 5.81e-02 7.91e-03 8.51e-01
5 1.8613e+00 1.8904e+00 1.57e-02 1.63e-03 1.58e-03 3.50e-02 3.10e-03 6.74e-01
6 2.1699e+00 2.1913e+00 9.85e-03 5.93e-04 4.34e-04 2.40e-02 8.88e-04 8.03e-01
7 2.5021e+00 2.5220e+00 7.93e-03 2.14e-04 1.12e-04 2.11e-02 2.37e-04 8.01e-01
8 2.7896e+00 2.8141e+00 8.76e-03 1.35e-04 4.21e-05 2.54e-02 9.20e-05 8.08e-01
9 3.2619e+00 3.2812e+00 5.91e-03 4.58e-05 9.63e-06 1.96e-02 2.16e-05 7.88e-01
10 3.6591e+00 3.6805e+00 5.86e-03 3.22e-05 3.07e-06 2.17e-02 7.16e-06 9.60e-01
11 4.0143e+00 4.0278e+00 3.36e-03 1.51e-05 1.27e-06 1.36e-02 3.07e-06 7.92e-01
12 4.3776e+00 4.3897e+00 2.76e-03 8.83e-06 5.17e-07 1.22e-02 1.31e-06 6.52e-01
13 4.8742e+00 4.8822e+00 1.65e-03 2.98e-06 1.20e-07 8.09e-03 3.36e-07 8.27e-01
14 5.0137e+00 5.0223e+00 1.72e-03 2.50e-06 7.68e-08 8.64e-03 2.22e-07 5.89e-01
15 5.3249e+00 5.3276e+00 5.03e-04 7.30e-07 2.35e-08 2.69e-03 6.77e-08 9.90e-01
16 5.4822e+00 5.4825e+00 6.11e-05 9.60e-08 3.22e-09 3.36e-04 9.26e-09 8.72e-01
17 5.5045e+00 5.5045e+00 1.14e-05 1.27e-08 4.29e-10 6.28e-05 1.23e-09 9.57e-01
18 5.5078e+00 5.5078e+00 6.45e-07 7.15e-10 2.43e-11 3.56e-06 6.98e-11 9.44e-01
19 5.5080e+00 5.5080e+00 5.11e-08 5.06e-11 1.74e-12 2.82e-07 4.99e-12 9.84e-01
20 5.5080e+00 5.5080e+00 1.21e-09 2.06e-12 5.94e-14 6.66e-09 1.18e-13 9.76e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 37.5ms
solution_summary(; result = 1, verbose = false)
├ solver_name : Clarabel
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 1
│ └ raw_status : SOLVED
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ ├ objective_value : -5.50801e+00
│ └ dual_objective_value : -5.50801e+00
└ Work counters
├ solve_time (sec) : 3.75230e-02
└ barrier_iterations : 20...and proves it by exhibiting the minimizer.
ν7 = moment_matrix(model7[:c])
η = atomic_measure(ν7, 1e-3)Atomic measure on the variables x[1], x[2] with 1 atoms:
at [2.3295201456728227, 3.1784931608225406] with weight 0.9999999549750614The first-order moments now coincide with the minimizer so round_solution finds it too:
round_solution(ν7, K, p)2-element Vector{Float64}:
2.3295200567388203
3.1784919250460106We can indeed verify that the objective value at x_opt is equal to the lower bound.
x_opt = η.atoms[1].center
p(x_opt)-5.508013306495363We can see visualize the solution as follows:
scatter!([x_opt[1]], [x_opt[2]], markershape = :star, label = nothing)
This page was generated using Literate.jl.